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Abstract 

The  paper  presents  calculations  for  the  radiation  of  a  loop  antenna  wrapped  around  a  metallic 
cylinder  into  a  conductive  medium.  The  problem  has  been  solved  analytically  by  using  Fourier 
transform  and  Bessel  functions.  The  evaluation  of  the  analytical  solution,  involves  the  use 
of  exponentially  scaled  Bessel  functions,  and  FFT  to  evaluate  inverse  Fourier  transform,  with 
all  computations  done  in  MATLAB.  The  analytical  solution  has  been  used  to  validate  a  goal- 
oriented  hp  adaptive  strategy  based  on  a  simultaneous  solution  of  a  dual  problem,  and  an  hp 
strategy  based  on  minimizing  the  projection-based  interpolation  error  of  a  reference  solution. 
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1  Introduction 

The  automatic  goal  driven  hp  adaptivity  has  been  used  to  solve  challenging  problems  in  computa¬ 
tional  electromagnetics  [2].  In  this  paper  we  formulate  and  solve,  both  analytically  and  numerically, 
a  problem  of  radiation  of  a  loop  antenna  wrapped  around  a  metallic  mandrell  into  a  conductive 
medium.  The  problem  has  been  solved  numerically  by  using  goal-oriented  hp  adaptive  code  [2] 
with  a  solution  evaluated  at  a  sequence  of  points.  The  main  goal  of  this  paper  is  to  solve  the 
same  problem  analytically  in  terms  of  Bessel  functions  and  inverse  Fourier  transform,  to  verify  the 
goal-oriented  hp  adaptive  computations.  Main  challenges  include: 

•  The  dynamic  range  of  the  solution  (ratio  between  the  maximum  and  minimum  values  of  the 
solution  within  the  domain  of  interest)  spans  14  digits. 

1On  leave  from  AGH  University  of  Science  and  Technology,  Department  of  Computer  Methods  in  Metallurgy, 
Cracow,  Poland ,  e-mail:  maciek@ices.utexas.edu,  mpaszynski@hotmail.com 
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•  In  order  to  evaluate  the  analytical  solution,  one  must  overcome  difficulties  with  the  numerical 
evaluation  of  the  (inverse)  Fourier  Transform  by  using  FFT.  This  requires  the  use  of  a  smooth 
truncating  function  to  enable  exponential  convergence  of  the  FFT  based  quadrature. 

The  presented  considerations  are  relevant  to  modeling  of  logging  devices  [3] .  The  high  accuracy  of 
solutions  for  this  class  of  problems  is  required  in  the  area  of  logging  while  drilling  inverse  problems. 
The  presentation  is  done  in  a  tutorial  style,  explaining  all  mathematical  details,  with  beginning 
graduate  students  in  mind. 


2  Statement  of  the  problem 


The  following  radiation  problem  is  relevant  to  drilling  technologies  [3].  A  loop  antenna,  wrapped 
around  an  infinite  metallic  cylinder,  radiates  into  a  conductive  homogeneous  medium,  as  illustrated 
in  Fig.  1.  The  radius  of  the  metallic  mandrell  is  a,  and  there  is  a  non-zero  distance  b  —  a  between 
the  metallic  mandrell  and  the  loop  antenna. 

We  start  with  the  time- harmonic  Maxwell  equations, 


V  x  Qv  x  E  +e 


u  n  yx,y,z,i) 


dt 2 


+  cr- 


dt 


dt 


where  y  stands  for  a  permeability,  cr  is  conductivity,  e  denotes  permittivity,  and  J'imp  is  an  imposed 
current. 


Fourier  transform  with  respect  to  time  t  leads  to, 


V  x  V  x  E  ( x ,  y,  z,  a;)^  —  (ui2e  —  iujo^j  E  ( x ,  y,  z,  u>)  =  —iuJimp  ( x ,  y,  z,  ix) .  (2-2) 


We  switch  now  to  the  axisynmietric  formulation  in  the  cylindrical  system  of  coordinates  with  z 
axis  located  on  the  axis  of  the  metallic  mandrell,  and  the  origin  at  the  center  of  the  loop  antenna. 
The  only  non-zero  component  of  the  electric  field  is  the  angular  component  E^.  Components 
Er  =  Ez  =  0,  due  to  the  axisymmetry  of  the  problem.  Maxwell  equations  (2.2)  reduce  to  a  single 
elliptic  problem, 


1  (  d2E(f>  1  _  IdE* 

y  y  dz 2  r 2  ^  r  dr 


—  (iu2e  —  aiujJ 


E^ 


(2.3) 


We  assume  that  the  imposed  current  is  given  by  a  Dirac  delta  Jimp  (. x ,  <fi,  z)  =  5  (r  —  b)  5  (z  —  0)  e^, 
where  is  the  angular  unit  vector  in  the  cylindrical  system  of  coordinates. 

For  the  numerical  computations  we  have  set  a  =  0.0254  [nr],  b  =  0.03048  [nr],  u>  =  2irf  where 
the  frequency  /  =  2- 106  [Hz],  y  =  yo  =  47rl0-7  is  the  free  space  permeability  and  e  =  eo  =  ^ylO^9 
is  the  free  space  permittivity.  The  metallic  mandrell  is  assumed  to  be  a  perfect  conductor.  The 
conductivity  of  the  medium  is  a  =  1  [S/m]. 
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The  formulation  can  be  transformed  into  a  non-dimensional  form,  by  defining, 

H  e  x  E 

—  €r  —  %r  —  ~j  Ej r  —  ,  (2-4) 

IMi  eo  d  Eq 

where  fir  is  the  dimensionless  permeability,  er  is  the  dimensionless  permittivity,  d  =  b  —  a  is  the 
spatial  scale  and  Eo  =  1  is  the  solution  scale.  These  are  the  four  basic  scaling  factors.  The 
remaining  quantities  are  transformed  into  a  dimensionless  form  by  defining, 


n  =  uiy/eoflod 


jimp 


jimp 


Eo 


(2.5) 


where  O  is  a  non-dimensional  wave  number,  S  is  a  non-dimensional  conductivity,  and  is  a 
non-dimensional  imposed  current. 

A  potential  pitfall  is  to  try  to  solve  the  equation  (2.3)  by  using  a  separation  of  variables  tech¬ 
nique.  There  are  two  issues: 


•  The  separation  of  variables  technique  can  be  applied  in  a  rigorous  way  only  to  a  problem 
resulting  in  a  Sturm-Liouville  system  [4] .  This  requires  the  domain  to  be  bounded  in  at  least 
one  direction.  Our  domain  of  consideration  is  unbounded  in  both  r  and  z. 

•  Equation  (2.3)  is  a  distributional  equation,  with  a  regular  distribution  part  associated  with 
the  left-hand  side,  and  a  non-regular  Dirac  delta  distribution  on  the  right-hand  side. 


We  begin  by  applying  the  Fourier  transform  with  respect  to  z  variable. 

The  Fourier  transform  of  a  non-regular  distribution  is  defined  by  generalizing  the  classical 
Fourier  transform.  We  start  from  the  Plancherel  theorem  [11]  which  states  that  the  L 2  scalar 
product  of  two  functions  is  equal  to  the  L 2  scalar  product  of  its  Fourier  transforms.  This  can  be 
interpreted  in  terms  of  distributions  as, 


<  FT  (Rf) ,  4>  >=<  Rf,  FT~l  (</>)  >, 


(2.6) 


where  FT  (Rf)  is  the  Fourier  transform  of  a  regular  distribution  Rf  associated  with  function  /,  0 
is  the  Fourier  transform  of  testing  function  (/>,  and  FT  V  (^j  is  the  inverse  Fourier  transform  of  (/>. 
Applying  (2.6)  to  the  Dirac  delta,  we  get, 


<  FT(5)A>=<  6,  FT~1  >=  FT"1  (0)  j  =  --L  J 4>  (u)  du  =<  R_l=J>, 

with  R_i_  denoting  the  regular  distribution  associated  with  factor  l—. 

s/2n  V27T 

Fourier  transforming  with  respect  to  z  variable,  we  obtain, 


E<f,  (r,  z)  - a 2  +  -2 - 

'  V  fY*  Zj  /y» 


1  dEf,  (r,  z)  8 2E4>  (r,  z) 


dr 


dr 2 


=  -iiiw-j=8  (■ r-b ) , 
V  27T 


(2.7) 


(2.8) 
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where  we  have  denoted  —a2  =  z2  —  \i  (J2 e  —  crico) .  All  terms  in  (2.8)  must  be  understood  as 
regular  distributions,  except  for  the  Dirac  delta  term,  and  all  derivatives  must  be  understood  in 
the  distributional  sense. 

With  the  Dirac  delta  concentrated  at  point  b  in  the  domain,  the  distributional  derivative  of  a 
regular  distribution  is  a  sum  of  a  regular  distribution  associated  with  the  classic  derivative 
E':  =/?£,,,  and  a  jump  term  [E^  (6)]  times  a  Dirac  delta  at  point  b  [6], 

V  -C/0 

-^  =  R^  +  [E<t>(b)]6(r~b).  (2.9) 

Consequently,  the  second  distributional  derivate  equals, 

=  +  (2-10) 

The  first  term  denotes  now  the  distributional  derivative  of  a  regular  distribution  associated  with  E ^ 
and  the  second  term  equals  the  jump  times  a  dipole  [6].  Absence  of  a  dipole  term  on  the  right-hand 
side  of  (2.8)  indicates  that  the  second  term  vanishes,  i.e.  E^  is  continuous  at  b.  Equation  (2.8) 
splits  now  into  tow  equations: 


•  the  regular  distribution  part, 

A  (-«2  +  ^)  -  -  E'i  =  0; 


•  the  singular  distribution  part 

—5  (r  —  b)  [E'^b)] 


b). 


(2.11) 


(2.12) 


Equation  (2.11)  is  the  Bessel  equation  for  which  a  solution  argument  must  be  scaled  by  factor  a 
[7].  We  split  our  ID  domain  into  two  intervals  (a,  b)  and  (6, +oo).  The  solution  over  each  interval 
is  a  linear  combination  of  the  Hankel  functions  of  the  first  and  the  second  type, 

E^b)  =  (or)  +  C2H[2)  ( ar ) ,  (2.13) 

E^+oo)  =  C3h[1]  (ar)  +  C4h[2)  (ar) .  (2.14) 

We  need  to  determine  now  the  four  coefficients  C\ ,  C2,  C3  and  C\. 

The  first  condition  to  be  used  is  the  jump  condition  (2.12). 

Moreover,  we  assume  that  the  metallic  mandrell  is  a  perfect  conductor.  With  the  electric 
held  vanishing  in  a  perfect  conductor,  and  the  tangential  component  of  E  being  continuous  across 
material  interfaces,  the  tangential  component  of  E  on  a  boundary  adjacent  to  a  perfect  conductor 
must  vanish  [8]  n  X  E^  =  0.  In  our  case  this  simply  translates  into  the  condition, 

E{^b)  (a)  =  0.  (2.15) 
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(2.16) 


The  third  condition  enforces  the  continuity  of  E ^  at  b, 

E^b)  ( b )  =  E^+oo)  ( b ) . 

The  fourth  condition  is  associated  with  the  asymptotic  behavior  of  Hankel  functions  [5].  It 
follows  that  the  Hankel  function  of  the  first  order  h\1'1  (a  (z)  r )  tends  to  0  when  the  imaginary  part 
of  a(z)  goes  to  infinity,  whereas  the  Hankel  functions  of  the  second  order  H\  ’  (a  (z)  r)  blows  up 
to  infinity,  when  the  imaginary  part  of  a  (z)  goes  to  infinity.  This  implies  that  C4  =  0. 

3  Closed  form  solution 


We  solve  now  the  system  of  three  equations  (2.12),  (2.15)  and  (2.16)  with  three  unknowns  C 1,  C2, 
C3.  The  solution  is, 


Cl  =  -C2 


H (aa) 
H (aa) 


(  (2) 

°3  ~  ~TJT  l"1 


(gj2)  (at)  -  )°°V 11  (at))  . 

\  H\  (aa)  ) 


Therefore,  the  angular  component  of  the  electromagnetic  field  in  the  Fourier  domain  is, 
M)  ,  ,  _  /  (1)  (2)  .  ,  H[1}  (ar)  h[2)  (aa)  h|1}  (ab)  \ 

E *  vrr  {ab)Hl  {ar) - 1^) - J’ 


(3.17) 


(3.18) 


(3.19) 


(3.20) 


£f+“)  M  =  (fl,(2)  (ab)  »!'’  ( ar )  -  H\  ^a) H «  (ab)  flj‘>  (ar) 

\  H{  (aa) 


(3.21) 


In  deriving  (3.18),  (3.19)  and  (3.20)  we  have  used  the  Wronskian  formula  for  Hankel  functions  [10], 


Hi  (ar)  H2  (ar)  —  H[  (ar)  H2  (ar)  = - . 

We  use  the  scaled  Hankel  functions,  in  order  to  eliminate  indeterminate  symbols, 

h[^  (z)  =  H j1)  (z)  exp  (—zi) , 

(z)  =  H[ 2'>  (z)  exp  (zi) . 

The  final  formulae  for  the  solution  in  the  Fourier  domain  are, 


ETb)  (r)  =  - 


(h[^  (ab)  H ^  (ar)  exp  (— ai(r  —  b )) 
4y  2  ' 


(3.22) 


(3.23) 

(3.24) 
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2 

y/ri 

2 

y/^2 

2 

— oo 

Too 

7T 

2 

Too 

7T 

2 

0 

y/ro 

00 

2 

y/ro 

00 

2 

Too 

Too 

7T 

2 

Too 

7T 

2 

Table  1:  Two  branches  of  a  (z)  =  +/  —  ^/rexp?}}. 


ifj  (aa) 

where  a  <  r  <  b,  and, 


where  a  <  b  <  r. 


^(b,+oo)  ^  ^14}2)  (ab)  H j1'*  (ar)  exp  (— m(r  —  6)) 

fr( 2)  /  \  \ 

—  Q°  (afr)  If}1)  (ar)  exp  (— a?(6  +  r  —  2a))  J  , 
if}  (aa)  J 


(3.25) 


(3.26) 


4  Numerical  implementation 

We  have  used  the  goal-oriented  hp  adaptive  code  to  compute  the  solution  at  a  sequence  of  points 
(r  =  b,z  =  zn),  where  zn  €  (0.5,  3.5)  meters.  The  point  source  has  been  replaced  with  a  small  finite 
size  antenna.  As  the  evaluation  of  the  analytical  solution  becomes  unstable  for  r  close  to  b,  we 
have  replaced  r  =  b  with  r  =  a  +  (b  —  a)  *  0.99.  Fig.  2  presents  the  real  and  imaginary  parts  of  the 
solution  as  a  function  of  I  for  the  selected  value  of  r. 

There  are  two  essential  components  of  the  numerical  evaluation  of  the  analytical  solution: 

•  numerical  evaluation  of  the  scaled  Hankel  functions, 

•  numerical  evaluation  of  the  inverse  Fourier  transform. 

Both  steps  have  been  implemented  in  MATLAB. 

4.1  Evaluation  of  Bessel  functions 

The  argument  for  all  Bessel  functions  of  the  third  kind  (Hankel  functions)  assumes  complex  values 
ar,  where  a  =  ( k 2  —  z2) 5  where  k 2  =  p,uj2e  —  i-iaui.  There  are  two  branches  of  the  complex  root 
+/  -  \J\k2  —  z2 1  exp?'}},  where  9  depends  upon  z,  as  illustrated  in  Fig.  3  and  Table  1.  The  selection 
of  the  branch  affects  the  asymptotic  behavior  of  the  solution  in  the  Fourier  domain  as  \z\  — >  oo,  and 
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(<2\  ^ 

it  corresponds  to  the  elimination  of  ;  (a  (z) see  the  discussion  in  the  previous  section.  The 
Hankel  functions  are  computed  by  using  MATLAB  routine  besselh(fi  ,K  ,Z  ,1  )  [9]  using  algorithms 
described  in  [10]. 


4.2  Evaluation  of  Fourier  transform  using  FFT 


The  inverse  FFT  is  used  to  approximate  the  integral, 


1  r+°o 

E<t>  ^  ^  =  2vr  J  X  ^  E<t>  ('r’  ^  eXP  ^ iz ^ 


(4.27) 


at  a  sequence  of  points,  where  x  denotes  a  truncating  function  defined  below.  The  integral  is 


computed  numerically  by  selecting  a  window 
approximation, 


-6,6 


in  the  Fourier  domain,  and  using  a  discrete 


2  N 


1  r+°  -  i 

E<t>  (r,  z)  ~  J  -  x  (z)  E'P  (r,  z)  exp  (izz)  dz  ^  —  2J  X  (%)  ^  (r>  %)  exp  (iz%)  dz, 
with  a  sequence  of  points, 


i=i 


Az 

zj  =  —  ~  N Az  +  ( j  -  1)  Az, 


for  j  =  1,  ...,2a.  The  inverse  FFT, 


2N 


(r,  Zk)  =  —J2x  ( zj )  (r,  %)  exp  (*zfc%)  dz, 


3= 1 


is  computed  at  a  sequence  of  points, 


Az 


Zfc  =  — - - NAz  +  (A;  —  1)  Az. 

The  size  of  the  real  domain  window  [—6,  6]  can  be  computed  from  relations, 

AzAz  =  —  and  66  =  nN. 

N 

There  are  three  general  observations  concerning  the  accuracy  of  the  inverse  FFT 


(4.28) 


(4.29) 


(4.30) 


(4.31) 


(4.32) 


A  twofold  increase  of  the  window  size  in  the  Fourier  domain 

L 

of  points  N ,  causes  the  real  domain  window  to  decrease  twofold  [— i6,  lb 


—26,  26j ,  under  a  fixed  number 
-6  -l 

2ui  2l 


A  twofold  increase  of  the  number  of  points  2 N,  under  a  fixed  size  of  the  Fourier  domain 
window,  causes  the  real  domain  window  to  increase  twofold  [—26,26], 

A  simultaneous  twofold  increase  of  the  window  size  in  the  Fourier  domain  —26,26  ,  with 
the  number  of  points  2 N  doubled,  increases  the  accuracy  with  the  real  domain  window  size 
remaining  fixed. 
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The  convergence  of  the  inverse  FFT  in  the  third  regime  is  illustrated  in  Fig.  4. 


Two  choices  of  the  truncating  function  x  have  been  investigated.  Results  obtained  with  a  simple 
characteristic  function  x  =  l[_g  ^  converge  very  slowly  and  do  not  match  the  FE  computations,  see 
Fig.  5. 

The  reason  is  that  the  integrated  function  along  with  its  derivatives  must  be  smooth  over  the 
complex  circle. 

Following  the  ideas  of  Bruno  at  al.  [1],  we  truncate  the  Fourier  integral  with  a  C q°  function 
(see  Fig.  6)  given  by, 


X(z)  =  { 


exp 

1 

exp 


i  (*~Vb+^) 


L  (*~^~8)) 


for  z  G 


—b,  —b  +  e 


for  5  G  (^—b  +  e,b  —  e)  , 


(4.33) 


for  5  G 


b  —  e,b 


where  e  =  Application  of  the  truncating  smooth  function  (4.33)  results  in  a  fast  (exponential  ?) 
convergence  of  the  inverse  FFT,  as  presented  in  Fig.  7. 

The  results  obtained  with  the  smooth  truncating  function  match  perfectly  the  numerical  results 
obtained  with  the  goal-oriented  hp  adaptive  strategy,  see  Fig.  8. 


5  Conclusions 

The  paper  presents  a  “struggle”  with  a  numerical  evaluation  of  the  inverse  Fourier  transform.  We 
have  shown  that  truncating  the  solution  in  the  Fourier  domain  with  a  Cq°  function,  guarantees  a 
fast  and  efficient  evaluation  of  Fourier  integral  using  FFT.  The  resulting  integrand  is  C°°  over  the 
complex  circle,  which  enables  a  fast  (exponential  ?)  convergence. 

The  presented  results  provide  an  additional  proof  for  the  unprecedented  accuracy  obtainable 
with  the  goal-oriented  hp  adaptivity. 

A  copy  of  the  MATLAB  code  can  be  requested  from  the  first  author. 


References 

[1]  O.  Bruno,  C.  Geuzaine,  F.  Reitich  “A  new  high-order  high-frequency  integral  equation  method 
for  the  solution  of  scattering  problems” ,  Proceedings  of  the  20th  Annual  Review  of  Progress  in 
Applied  Computational  Electromagnetics ,  Syracuse,  New  York,  (2004) 

[2]  D.  Pardo,  L.  Demkowicz,  C.  Torres- Verdin,  L.  Tabarovsky,  “A  Goal  Oriented  /ip-Adaptive 
Finite  Element  Strategy  with  Electromagnetic  Applications.  Part  I:  Electrostatics.”,  Submitted 
to:  International  Journal  for  Numerical  Methods  in  Engineering,  (2004) 


[3]  J.  R.  Lovell,  Finite  Element  Methods  in  Resistivity  Logging,  Delft  University  of  Technology, 
(1993),  p.38-45 

[4]  M.D.  Greenberg,  Foundations  of  Applied  Mathematics ,  Prentice-Hall  International,  (1978) 

[5]  N.  N.  Lebedev,  Special  Functions  and  Their  Applications  ,  Dover  Publications,  New  York, 
(1972) 

[6]  A.H.  Zemanian,  Distribution  Theory  and  Transform  Analysis:  An  Introduction  to  Generalized 
Functions,  With  Application,  Dover  Publications,  New  York,  (1965) 

[7]  J.  S.  Walker,  Fourier  Analysis,  Oxford  University  Press,  New  York,  (1988) 

[8]  L.  Demkowicz,  “Finite  Element  Methods  for  Maxwell  Equations”,  Encyclopedia  of  Computa¬ 
tional  Mechanics,  John  Wiley  &  Sons,  (2004) 

[9]  MATLAB  User  Guide 

[10]  M.  Abramowitz,  I.  A.  Stegun,  Handbook  of  Mathematical  Functions,  National  Bureau  of  Stan¬ 
dards,  Applied  Math.  Series  #55,  Dover  Publications,  (1965) 

[11]  J.  T.  Oden,  L.  Demkowicz,  Applied  Functional  Analysis,  CRC  Press,  1996 


9 


1 

loop  antenna 

- 

cT 

J  imP  impressed  current 

a _ \  r 

c 

metallic  mandrel 

Figure  1:  Configuration  of  the  metallic  mandrel  and  loop  antenna 
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Figure  2:  Analytic  solutions  for  r  =  (a  +  (b  —  a)  *  0.99).  Red  color  denotes  a  real  part,  blue  color 
denotes  an  imaginary  part. 


Figure  3:  Two  branches  of  az  =  +/  —  yVexpi^.  Vertical  axis  denotes  real  part  of  az,  horizontal 
axis  denotes  imaginary  part  of  az. 
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Figure  4:  Convergence  of  the  inverse  Fourier  transform  using  FFT  without  truncating  function 
X  =  lj_ggp  for  N  =  25’6’  -’11,12  and  b  =  3000* 2— 7’  6,— .  Vertical  axis  =  logio (Electromagnetic 
field);  Horizontal  axis  =  distance  [m].  Red  color  denotes  a  real  part,  blue  color  denotes  imaginary 
part  of  the  solution.  Black  stars  denote  real  part  of  the  goal-oriented  hp  adaptive  result,  pink  stars 
denote  its  imaginary  part. 
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Figure  5:  Comparison  of  the  inverse  Fourier  transform  using  FFT  without  truncating  func¬ 
tion  x  =  1  [_g  £]  with  the  results  of  goal-oriented  hp  adaptive  computations.  Vertical  axis  = 
logio (Electromagnetic  field);  Horizontal  axis  =  distance  [nr].  Red  color  denotes  a  real  part,  blue 
color  denotes  imaginary  part  of  the  solution.  Black  stars  denote  real  part  of  the  goal-oriented  hp 
adaptive  result,  pink  stars  denote  its  imaginary  part. 


Figure  6:  Truncating  function  \ 
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Figure  7:  Convergence  of  the  inverse  Fourier  transform  using  FFT  with  the  smooth  truncating 
function  x  given  by  (4.33),  for  N  =  25>6>-’11>12  an(j  £  —  3000  *  2“7,”6,_5,***,“1,0.  Vertical  axis  = 
logio (Electromagnetic  field);  Horizontal  axis  =  distance  [m].  Red  color  denotes  a  real  part,  blue 
color  denotes  imaginary  part  of  the  solution.  Black  stars  denote  real  part  of  the  goal-oriented  hp 
adaptive  result,  pink  stars  denote  its  imaginary  part. 
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Figure  8:  Comparison  of  the  inverse  Fourier  transform  using  FFT  with  truncating  function  \ 
given  by  (4.33)  with  the  results  of  goal-oriented  hp  adaptive  computations.  Vertical  axis  = 
logio (Electromagnetic  held);  Horizontal  axis  =  distance  [nr].  Red  color  denotes  a  real  part,  blue 
color  denotes  imaginary  part  of  the  solution.  Black  stars  denote  real  part  of  the  goal-oriented  hp 
adaptive  result,  pink  stars  denote  its  imaginary  part. 
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